Carderock  Division,  Navai  Surface  Warfare  Center 

Bethesda,  Maryland  20084-5000 


(0 

c 

g 

3 

O 

cd 

O 

o 

cd 

X 

LU 


"D 

2 

CO 

Q. 

E 

o 

O 

0) 

CO 

c 

o 

Q. 

CO 

O 

cr 

O) 

TD 

O 


5  < 

O 
t 


CO 


0) 


0) 


T3 

<1  :i 


CD 


O 
O  O 
« 
o 

CCl  LU 


NSWCCD-TR-96/006  March  1996 

Survivability,  Structures  and  Materials  Directorate 
Research  and  Development  Report 


Ref  iected-Afterf  iow  Virtual-Source  (RAVS) 
Model  Response  Compared  to  Exact 
Calculations  for  Elastic  Cylinders 
Attacked  by  Planar  Waves 


by 

George  V.  Waldo,  Jr. 


19960402  129 


Approved  for  public  release;  distribution  is  unlimited. 


UNCLASSIFIED 

SECURITY  CLASSiFICATIQN  OF  THIS  PAGE 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
0MB  No  0704  0188 


la  REPORT  security  CLASSIFICATION 

UNCLASSIFIED 


lb  RESTRICTIVE  MARKINGS 


2a  SECURITY  CLASSIFICATION  AUTHORITY 


2b.  DECLASSIFICATION /DOWNGRADING  SCHEDULE 


3  DISTRIBUTION/AVAILABILITY  OF  REPORT 

Approved  for  public  release;  distribution 
is  unlimited. 


4.  PERFORMING  ORGANIZATION  REPORT  NUMBER(S) 

NSWCCD-TR- 9 6 / 0 0 6 


5  MONITORING  ORGANIZATION  REPORT  NUMBER(S) 


6a.  NAME  OF  PERFORMING  ORGANIZATION 

Naval  Surface  Warfare  Center 
Carderock  Division 


6b  OFFICE  SYMBOL 
(If  applicable) 


7a.  NAME  OF  MONITORING  ORGANIZATION 


Code  672 


6c  ADDRESS  {City,  State,  and  ZIP  Code) 


7b  ADDRESS  (City,  State,  and  ZIP  Code) 


Bethesda,  MD  20084-5000 


8a.  NAME  OF  FUNDING /SPONSORING 
„  ORGANIZATION  „  , 

Naval  Surface  Warfare  Center 
Carderock  Division 


8b.  OFFICE  SYMBOL 
(If  applicable) 

Code  on 


9  PROCUREMENT  INSTRUMENT  IDENTIFICATION  NUMBER 


8c  ADDRESS  (C/fy,  State,  and  ZIP  Code) 


10  SOURCE  OF  FUNDING  NUMBERS 


PROGRAM 

PROJECT 

TASK 

WORK  UNIT 

ELEMENT  NO 

NO 

NO 

ACCESSION  NO 

1 1 ,  TITLE  (Include  Security  Classification) 

Ref lected-Af terf low  Virtual-Source  (RAVS)  Model  Response  Compared  to  Exact 
Calculations  for  Elastic  Cylinders  Attacked  by  Planar  Waves _ 


12.  PERSONAL  AUTHOR(S) 

Waldo.  Jr..  George  V, 


13a.  TYPE  OF  REPORT 
Final 


13b  TIME  COVERED 

FROM  950315  TO950715 


14.  DATE  OF  REPORT  {Year,  Month,  Day) 

96  February 


15  PAGE  COUNT 
60 


16  SUPPLEMENTARY  NOTATION 


1  17.  COSATI  CODES  | 

FIELD 

GROUP 

SUB-GROUP 

18  SUBJECT  TERMS  {Continue  on  reverse  if  necessary  and  identify  by  block  number) 

acoustic  virtual-source  shell 

reflection  elastic  cylindrical 

afterflow  structural  approximation 


19  ABSTRACT  {Continue  on  reverse  if  necessary  and  identify  by  block  number) 

A  recently  developed  theoretical  approximation,  the  Ref lected-Af terf low  Virtual- 
Source  (RAVS)  model,  is  applied  to  the  case  of  a  planar  shock  wave  attacking  an  elastic 
cylindrical  shell.  To  do  this,  an  equation  of  motion  for  a  structural  surface  is 
presented.  Also,  equations  for  the  pressure  that  develops  when  an  acoustic  wave 
interacts  with  a  curved  and  compliant  surface  are  presented.  These  equations  were 
derived  using  the  RAVS  model.  An  expression  for  the  velocity  of  the  structural 
surface  is  derived  from  the  equation  of  motion  and  a  numerical  method  is  developed 
using  this  expression.  This  is  applied  to  the  case  of  a  planar  wave  attacking 
an  elastic  cylindrical  shell.  Calculations,  using  the  RAVS  model,  are  shown  to  be 
in  good  agreement  with  all  of  the  exact  calculations  that  were  published  by  Huang 
for  all  locations  and  times .  It  is  concluded  that  the  RAVS  model  gives  a  good 
approximation  for  the  loading  and  response  of  an  elastic  cylindrical  shell  attacked 
by  a  planar  wave.  The  RAVS  model  is  shown  to  be  even  better  for  a  spherical  wave 


20.  DISTRIBUTION /AVAILABILITY  OF  ABSTRACT 
□  UNCLASSIFIED/UNLIMITED  El  SAME  AS  RPT 


□  DTIC  USERS 


21. 


ABSTRACT  SECURITY  CLASSIFICATION 

UNCLASSIFIED 


22a 


NAME  OF  RESPONSIBLE  INDIVIDUAL 

George  V.  Waldo,  Jr. 


22b  TELEPHONE  (/nc/ude  Area  Code) 

(301)  227-1782 _ 


22c  OFFICE  SYMBOL 

Code  672 _ 


_ UNCLASSIFIED _ 

SECURITY  CLASSIFICATION  OF  THIS  PAGE 


CONTENTS 

Page 

ABSTRACT . 1 

INTRODUCTION . 1 

EQUATION  OF  MOTION . 2 

PRESSURE  OF  THE  SCATTERED  WAVE . 2 

PRESSURE  OF  THE  RADIATED  WAVE . 5 

SOLUTION  OF  THE  EQUATION  OF  MOTION . 6 

NUMERICAL  CALCULATIONS . 7 

EXPONENTIALLY  DECAYING  WAVE . 9 

PLANAR  ATTACKING  WAVE . 11 

CYLINDRICAL  STRUCTURAL  SURFACE . 11 

COMPARISONS  TO  EXACT  CALCULATIONS . 13 

Situation  with  M=2,  h/a^=l/31  and  C^=  13.685665  . 14 

Situation  with  M=4.4189,  h/ai=l/69  and  0^=12.581197  . 16 

Situation  with  M=6.  41975,  h/ai=l/100,  0^=13 . 685665 , 

and  g=0 . 20 

Situation  with  M=9. 09184,  h/ai=l/142,  0^=12.581197, 

and  q=0 . 21 

Situation  with  M=4.4189,  h/ai=l/69,  0^=12.581197, 

and  q=4 . 8c /a^ . 22 

RIGID  CYLINDERS . 23 

SUMMARY  AND  CONCLUSIONS . 25 

APPENDIX  A.  SOLUTION  OF  THE  EQUATION  OF  MOTION . 27 

APPENDIX  B .  STRUCTURAL  RESISTANCE  PRESSURE . 33 

APPENDIX  C  .  COMPUTER  PROGRAM . 37 

REFERENCES .  . . 47 

iii 


FIGURES 


Page 


1.  Geometry  for  virtual  source . 4 

2  .  Planar  wave  attacking  cylindrical  shell . 12 

3.  Nondimens ional  velocity  vs.  time,  0=0,  g=0,  M=2 . 15 


g=0,  M=4.4189 
g=0,  M=4.4189 
q=0,  M=4.4189 
g=0,  M=4.4189 
g=0,  M=4.4189 
g=:0,  M=4.4189 


4.  Nondimens  ional  displacement,  0-th  mode,  q=0,  M=2 .  .  . 

5.  Nondimens  ional  velocity  vs.  time,  1-st  mode,  g=0,  M=2 .  .  . 

6.  Nondimens  ional  velocity  vs.  time,  0=0,  g=c/a^,  M=2 .  .  . 

7.  Nondimens  ional  velocity  vs.  time,  0=0,  g=0,  M=4.4189 

8.  Nondimensional  velocity  vs .  time,  0=71/2, 

9.  Nondimensional  velocity  vs.  time,  0=7C, 

10.  Strain  vs.  nondimensional  time,  0=0, 

11.  Strain  vs.  nondimensional  time,  0=7t/2, 

12.  Strain  vs.  nondimensional  time,  0=7i, 

13.  Nondimensional  displacement,  0-th  mode, 

14.  Nondimensional  velocity  vs .  time,  1-st  mode,  g=0,  M=4.4189 

15.  Strain  vs.  nondimensional  time,  0=0,  q=3c/a^,  M=4.4189 

16.  Strain  vs.  nondimensional  time,  0=71/2,  q=3c/a^,  M=4.4189 

17.  Strain  vs.  nondimensional  time,  0=7i,  q=3c/a^,  M=4.4189 

18.  Nondimensional  displacement,  0-th  mode,  M=6 .  41975 ...  . 

19.  Nondimensional  velocity  vs .  time,  1-st  mode,  M=6 . 41975 ...  . 

20.  Nondimensional  displacement,  0-th  mode,  M=9.  09184.. 

21.  Nondimensional  velocity  vs.  time,  1-st  mode,  M=9. 09184.. 

22.  Nondimensional  velocity  vs .  time,  0=0,  q=4.8c/a^,  M=4.4189..22 

23.  Nondimensional  velocity  vs .  time,  Q=n/6,q=4 . 8c/aj,  M=4.4189..22 

24.  Nondimensional  velocity  vs .  time,  0=%7C,  q=4.8c/ai,  M=4. 4189.. 23 

25.  Nondimensional  velocity  vs .  time,  0=%7C,  q=4.8c/ai,  M=4. 4189.. 23 

26.  Nondimensional  velocity  vs.  time,  0=7C,  q=4.8c/aj^,  M=4. 4189.. 23 


.15 

.16 

.16 

.16 

.17 

.17 

.18 

.18 

.18 

.19 

.19 

.20 

.20 

.20 

.20 

.21 

.22 

.22 


XV 


27.  Rigid  cylinder  attacked  by  a  spherical  wave  with  constant 


pressure,  nearest  distance  =  0.5ai,  longitudinal  distance  =  0 . .  .24 

28.  Rigid  Cylinder  attacked  by  a  spherical  wave  with  constant 

pressure,  nearest  distance  =  0.5ai,  longitudinal  distance  =  1.41421ai . 24 


V 


THIS  PAGE  IS  INTENTIONALLY  BLANK. 


vi 


ABSTRACT 


A  recently  developed  theoretical 
approximation,  the  Ref lected-Afterf low 
Virtual-Source  (RAVS)  model,  is  applied  to 
the  case  of  a  planar  shock  wave  attacking  an 
elastic  cylindrical  shell.  To  do  this,  an 
equation  of  motion  for  a  structural  surface 
is  presented.  Also,  equations  for  the 
pressure  that  develops  when  an  acoustic  wave 
interacts  with  a  curved  and  compliant  surface 
are  presented.  These  equations  were  derived 
using  the  RAVS  model.  An  expression  for  the 
velocity  of  the  structural  surface  is  derived 
from  the  equation  of  motion  and  a  numerical 
method  is  developed  for  using  this 
expression.  This  is  applied  to  the  case  of  a 
planar  wave  attacking  an  elastic  cylindrical 
shell.  Calculations,  using  the  RAVS  model, 
are  shown  to  be  in  good  agreement  with  all  of 
the  exact  calculations  that  were  published  by 
Huang  for  all  locations  and  times.  It  is 
concluded  that  the  RAVS  model  gives  a  good 
approximation  for  the  loading  and  response  of 
an  elastic  cylindrical  shell  attacked  by  a 
planar  wave.  The  RAVS  model  is  shown  to  be 
even  better  for  a  spherical  wave  attack  which 
is  a  case  of  more  practical  interest. 


INTRODUCTION 

Calculation  of  the  deformation  of  a  structure  caused  by  the 
pressure  wave  from  an  underwater  explosion  is  of  much  interest  to 
the  Navy.  To  gain  this  capability  an  approximate  model  was 
introduced^' ^  to  determine  the  pressure  that  develops  when  an 
acoustic  wave  interacts  with  a  curved  and  compliant  surface. 

This  model  uses  the  virtual-source  concept  in  analogy  with 
geometrical  optics.  Because  this  model  includes  the  afterflow 
velocities  of  both  the  attacking  and  reflected  waves,  it  is 
called  the  "Ref  lected-Afterf  low  Virtual-Source"  (RAVS)  model.  In 
this  paper,  an  equation  of  motion  for  a  structural  surface  is 
obtained  using  the  RAVS  model.  From  this  equation  of  motion,  the 
deformation  of  the  surface  is  calculated.  A  computer  program  is 


1 


presented  to  determine  the  deformations  for  the  special  case  of  a 
planar  shock  wave  interacting  with  a  cylindrical  elastic  shell. 
Calculations  with  this  program  (using  RAVS)  are  shown  to  be  in 
good  agreement  with  the  exact  calculations  by  Huang. 

EQUATION  OF  MOTION 

The  equation  of  motion  for  the  velocity  u  (t)  of  the  surface 
in  the  direction  that  is  opposite  to  the  normal  to  the  surface 
(Figure  1)  is 

muit)  =p^^^(t)  -p^^j.(t)  ,  (1) 

where  m  is  the  total  mass  (including  the  attached  stiffeners)  of 
the  structural  surface  per  unit  area,  t  is  the  time  after  arrival 
of  the  attacking  wave,  p^a^(t)  is  the  total  pressure  in  the  fluid, 
and  Pstr(t)  is  the  structural  pressure,  i.e.,  the  force  per  unit 
area  due  to  the  structure  in  the  direction  that  is  normal  to  the 
structural  surface  (pointing  into  the  fluid) .  This  quantity  is 
determined  by  the  properties  of  the  structure.  The  total 
pressure  in  the  fluid  at  the  point  of  interest  on  the  surface  of 
the  structure  is 

Ptot  (  t)  =  Pftydro  P  ( t)  +  ( t)  +  >  (2 ) 

where  Phydro  is  the  hydrostatic  pressure,  p(t)  is  the  pressure  of 
the  attacking  wave,  is  the  pressure  of  the  scattered 

wave,  and  Pr^^(t)  is  the  pressure  of  the  radiated  wave. 


PRESSURE  OF  THE  SCATTERED  WAVE 

Because  the  scattered  wave  is  defined  to  be  the  wave  that 
would  occur  if  the  surface  did  not  move  (if  there  were  no 
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separation  of  the  fluid  from  the  surface)  ,  the  velocities  of  the 
attacking  and  scattered  waves  in  the  direction  of  the  normal  to 
the  surface  must  be  equal,  that  is, 


t 

^  (pit')dt' 

COS0  = 

pc  prJ^ 

Pscat(t) 


pc 


-pJPscat(t')  dt' 


COsP  , 


(3) 


which  implies  that  (see  Eq.  13  of  Ref.  1)  the  pressure  of  the 
scattered  wave  is 


Pscat  ^ 


pit)  *  /  ®  ^  pit')  dt' 


cos6 


COSp 


(4) 


where  0  (see  Figure  1)  is  the  angle  of  incidence  of  the  attacking 
wave  ( i .  e . ,  the  angle  between  the  normal  to  the  surface  and  the 
line  from  the  source)  ,  r  is  the  distance  from  the  source,  p  is 
the  mass  density  of  the  fluid,  c  is  the  speed  of  sound. 


b  =  ^b^^  *  Aaia-b^)  [sin(o/2)]=^ 


(5) 


is  the  distance  to  the  virtual  source,  and  the  cosine  of  the 
angle  of  the  scattered  wave  P  (between  the  radial  direction  of 
the  equivalent  sphere  and  the  line  to  the  virtual  source)  is 


bo  +  2  (a -bo)  [sin  (a/2)  ]  2  . 

COSp  =  — - ^ ' 

b 

In  these  equations,  the  radius  of  the  equivalent  sphere  at  the 
point  of  interest,  a,  is  given  by 


1 

a 


l{jL 

2(a, 


(7) 


where  a^  is  the  radius  of  curvature  in  a  given  direction  on  the 
actual  surface  and  a^  is  the  radius  of  curvature  in  the  direction 
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that  is  perpendicular  to  the  given  direction.  As  easily  shown, 
the  radius  of  this  equivalent  sphere  does  not  vary  on  the  choice 
of  the  direction  a^.  If  the  actual  surface  were  a  sphere,  then 
a^=a2=a.  However,  if  the  actual  surface  were  a  cylinder,  then  a^ 
is  infinite  and  a=2aj.  Also,  the  distance  bg  from  the  surface  of 
the  equivalent  sphere  to  this  virtual  source  is  given  by 


=  -i  +  —  ,  for  d  ^  0  and  a  >  0 , 
jDq  d  a 


(8) 


where 


d  = 


+  2Tacosd 


+  a^  +  2racos6  +  a 


(9) 


Figure  1.  Geometry  for  virtual  source. 
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is  the  distance  from  the  source  to  the  surface  of  the  equivalent 
sphere  .*  The  angle  between  the  line  from  the  source  to  the 
center  of  the  equivalent  sphere  and  the  line  from  the  center  of 
the  sphere  to  the  point  of  interest  on  the  structural  surface  is 


a  =  C  ,  for  a  +  rcos6  ^  0  , 


and 


a  =  C  +  71  /  for  a  +  rcos0  <  0  , 


where 


C  =  atanl^Si^] 
\  a  +  rcosB 


(10) 


(11) 


(12) 


PRESSURE  OF  THE  RADIATED  WAVE 

The  pressure  at  the  point  of  interest  due  to  the  radiated 
wave  also  is  modeled  as  if  it  were  emanating  from  the  virtual 
source.  This  pressure  is  determined  by  solving 


u(t)  =  - 


Prad(t) 


pc 


1 

p2> 


w 

-r/c 


cosP  , 


(13) 


where  u(t)  is  the  velocity  of  the  surface  in  the  direction  that 
is  opposite  to  the  normal  to  the  surface,  i.e.,  the  velocity  in 


*For  cases  with  d  <  0,  the  virtual  source  is  assumed  to  be 
located  at  the  actual  source.  This  implies  that 

b  =  r  and  P  =  ti  -  6  , 

see  the  discussion  after  Eq.  14  of  Ref.  1. 

Under  certain  conditions,  a  calculation  might  be  acceptable 
for  a  concave  surface  (i.e.,  a  surface  with  an  equivalent  sphere 
having  a  negative  radius  (7)),  However,  there  are  complications 
that  could  occur  for  other  conditions.  Thus,  for  the  models  that 
are  presented  in  this  report,  it  is  required  that  a  >  0.  For 
numerical  calculations,  if  a  were  negative,  then  it  would  have  to 
be  set  equal  to  a  very  large  positive  nvimber  so  that  the  effective 
curvature  in  the  calculation  would  be  approximately  zero. 
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the  direction  into  the  surface  of  the  structure.  (Note  that  the 
opposite  sign  convention  is  used  in  Huang,*)  Thus,  the  pressure 
of  the  radiated  wave  is 


p„s(t)  =-«■<■/  e  ^  ,  (1*) 

-z/c  " 

as  in  Eq.  22  of  Ref.  1. 

Of  course,  it  is  assumed  that  the  deformations  of  the 
surface  are  so  small  that  they  have  a  negligible  effect  on  the 
curvatures.  If  these  curvatures  were  continually  updated  during 
a  calculation  to  give  the  current  curvature,  there  should  be  an 
improvement  in  the  calculation. 


SOLUTION  OF  THE  EQUATION  OF  MOTION 

In  Appendix  A,  it  is  shown  that  the  solution  of  (1)  is 

u(t)  =  Uo(t)  -  Ugf.^{t)  ,  (15) 

in  which  the  velocity,  if  there  were  no  structural  pressure,  is 


Uoit) 


lit)  L  ^  cose\  ^  _!/ c 
m  \  cosP/  m\b 


+ 


c  cos9 
r  cosP 


lit) 


(16) 


and  the  retarding  velocity  due  to  the  structural  pressure  is 


m 


-f/c 


where 


lit)  =  e-’itj  e’it'p(tO  dt'  , 


(17) 


(18) 


6 


(19) 


and 


TncosP 


(20) 


is  the  damping  coefficient,  see  also  (82)  .  If  Pstr(^)  is  not 
known,  then  it  must  be  determined  from  the  deformation  of  the 
structural  surface. 

NUMERICAL  CALCULATIONS 


write 


To  calculate  I(t)  nxomerically  (see  (18)),  it  is  useful  to 

t+At 


J(  t  +  At) 


g-i|{t  +  At) 

g-t,Atj(t) 


(21) 


Thus,  the  increment  of  I(t)  is 


where 


AJ  =  J(t+At)  -  J(t) 


t+At 

=  g-H(t+At)  /  dt' 

t 


(l-e-’>A‘=)  j(t) 


oH(t+At)  _  o’lt 

«  — - ^p(t+At)  -  (1 -e‘^^‘')r(t) 

=  tl'Ml -e-^^*^)p(t+At)  -  (1 -e-"^^)  J(t) 

=  [p(  t+At)  -■nJ(t)  ]  Q  , 

(22) 


®  =  |exp(-|At)sinh(f  At) 


(23) 


This  would  be  exact  if  p(t)  were  constant  in  the  interval  from  t 
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to  (and  including)  t+At.  Note  that 


J(t)  =0  ,  for  t  i  0  . 

The  numerical  calculation  of  Iscr(^)  (see  (19))  can 
performed  in  a  similar  fashion.  Thus,  the  increment  in 
is 

^fatz  *  ^Patz  (t  +  At)  ~P]iydzo  ~  ^  ^atz  (  t)  ]  0 

and 


J^,^(t  +  At)  =  J^,^(t)  +AJ^,^  , 


where 


t^-r/c  . 


Differentiation  of  (17),  u^t.^(t),  gives 

Thus,  the  increment  of  Uscr(t)  is 


AU^tr  »  [Alatr*  latrif^'i  ic/b)  At]  /Ol 


^str  '^Btz 


and 


U^^^{t  +  At)  =  U^tr(t)  ' 

where 

Ugf.j.U)  -  0  ,  for  ti-r/c  . 

Thus,  the  velocity  of  the  surface  (that  is  in  the 
opposite  to  the  normal  to  the  surface)  is 

U(t)  =  Uo(t)  -Ugtj.{t)  , 

as  given  by  (15) . 


(24) 

be 

I.tr  (t) 

(25) 

(26) 

(27) 

(28) 

(29) 

(30) 

(31) 

direction 

(32) 
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EXPONENTIALLY  DECAYING  WAVE 

Now  consider  an  attacking  pressure  wave  with  a  sudden  rise 
and  an  exponential  decay, 


p ( t)  =  0  ,  for  t^O, 

=  t>  0  , 


(33) 


where  Po  is  the  initial  pressure  and  g  is  the  exponential  time- 
decay  constant.  For  this  case,  it  is  convenient  to  define 


s  =  iUlS 
2 


(34) 


and 


=  JLl3 


with  these  definitions  (18)  becomes 


I  t  =  Po- - - -  =  Poe  — 

°  rj  -  g  °  2x 

-  Pote'"*^8h(x^)  /  t>0  , 


(35) 


where 


sh{x)  =  ,  for  x*0  , 


(36) 


(37) 


=  1  ,  for  x=0  . 

Thus,  the  velocity  if  there  were  no  structural  pressure,  (16),  is 

U,(t)  =  =  w).  ,  (38) 

°  jn  \  cosp/  \b  r  cosP/  m 


where,  for  nximerical  calculations. 
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J(t)  =  ^fp(t')  dt'- 


n  Ti  -g 


_  Po  n  -  Q-r\e~'^  +  Qe'''^  _  ^  -  (s  +  x)  +  (s-x)  ‘ 

n  g{ti-g)  tigi  2% 

=  jE2.fl  -  c-se  (s  +  x)e^’=-  (s-x)e'^ 

^<71  2x 


=  ^li 


iig 

ng 


|l -e-"^  cosh(xt)  H-sSjrSjllXfc)  j 


{l  -  e'®^  [cosh(xt)  +  St  3h(xt)  ]}  ,  for  qt>  e.  and  i]  t  >  e  , 


(39) 


in  which 

1>>€>0,  (e.gr.  ,e  =  0.0001)  ,  (40) 

see  (37)  for  sh(x).  Otherwise,  for  numerical  calculations, 

1  -  e  1  -  e  -»i  t 


Jit)  =Po - 2. 


r\-g 


Po 


Ti-gL<? 


-  qt 

—  e  2  sinhi 


qt 
^2 


I  >1 


e  ^  sinh|-^| 


Pot 


-g| 


e-^s^f)-e 


(41) 


for  gt  ^  e  and/or  t)  t  5  e  and  >  (es)  ^  , 


otherwise,  for  numerical  calculations. 


•^(t)  =  Pq—  /  for  gt^e  and/or  T]t^e  and  ^  .  (42) 

/j 
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PLANAR  ATTACKING  WAVE 


For  a  planar  attacking  wave,  the  source  is  at  an  infinite 
distance, 


r  -  00  . 


In  this  case. 


d 


and 


a  -  6  / 
bo  a/2  , 


see  (8) -  (12)  . 


(43) 

(44) 

(45) 

(46) 


CYLINDRICAL  STRUCTURAL  SURFACE 

For  a  cylindrical  structural  surface,  32  is  infinite.  Thus, 
for  a  planar  wave  propagating  perpendicular  to  the  axes  of  the 
cylinder, 

a  =  2a^  and  bo  =  a^  ,  (47) 

see  (7)  and  (8)  .  Also,  the  angle  of  incidence,  0,  is  the  same  as 
the  angle  of  the  point  of  interest  on  the  cylinder  above  the 
direction  of  propagation  of  the  wave,  see  Figure  2. 

For  a  cylindrical  structural  surface  (i.e.,  shell.  Figure  2) 
composed  of  elastic  material  with  Young's  modulus,  E,  Poisson's 
ratio,  V,  thickness,  2h,  and  radius  a^,  it  is  shown  in  Appendix  B 
(see  (97))  that  the  structural  resistance  pressure  is 


P 


2  Eh 


str 


(l-v^)ai 


^radial 


J/L 


j(4) 

^zadial 


- 


(1) 


(3) 


(48) 
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where  d^adiai  is  the  deflection 
in  the  radial  direction,  de  is 
the  deflection  in  the 
tangential  direction,  and  the 
superscripts  (in  parentheses) 
indicate  partial 
differentiation  with  respect 
to  6.  Also,  it  is  shown  in 
Appendix  B  (see  (98))  that  the 
equation  of  motion  for  the 
tangential  deflection  is 


Figure  2.  Planar  wave  attacking 
cylindrical  shell. 


2  Eh 


m 


(l-v2)ai^ 


(‘•ii') 


- 


j(l) 

^radial 


radial 


(49) 


in  which  the  mass  per  unit  area  of  the  shell  is 

m  =  2pgh  ,  (50) 

where  is  the  mass  density  of  the  shell  material. 

For  numerical  calculations,  the  radial  deflection  is 
approximated  as 

(  t)  +  [u  ( t)  +  U  ( t +A  t)  ]  A  t/2  ,  (51) 

using  (32)  for  u(t)  and  (48)  for  Pst:r(^)  •  Using  (49),  the 
tangential  velocity  is  approximated  as 

d^(t  +  At)  »  <^(t)  +  c^(t)  At  (52) 

and  the  tangential  deflection  is 

c^{t  +  At)  »  c^(t)  +  [4(t)  +(^(t  +  At)]  At/2  .  (53) 
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A  computer  program  (written  in  the  Mathematical  language) 


using  these  equations  is  presented  in  Appendix  C. 


COMPARISONS  TO  EXACT  CALCULATIONS 

Exact  calculations  for  this  problem  were  performed  by  Huang^ 
with  Pt^y^^=0.  In  Eq.  (1)  of  Huang's  paper,  the  following 
dimensionless  quantities  are  introduced: 

M  =  pa^^/ {2hpg)  and  =  E/ [pg(i-\^)  c^]  .  (54) 

His  results  for  these  quantities  can  be  regarded  as  functions  of 
the  number  of  radial  transit  times  after  arrival  of  the  front  of 
the  attacking  wave  at  0=0,  i.e., 

T*  =  cT/a^  ,  (55) 

where  T  is  the  time  after  arrival  of  the  attacking  wave  at*  0=0. 

In  a  similar  fashion,  the  nondimens ional  radial  velocity  is 
defined  as 


n* 


(56) 


where 


Po 


(57) 


pc 

is  the  velocity  of  the  fluid  immediately  behind  the  front  of  the 
attacking  wave  (see  Cole,®  p.  35)  . 


'Thus,  the  time  after  arrival  of  the  front  at  a  given  location 
is 

t  =  r  -  —  (1 -COS0)  . 
c 
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Situation  with  M=2,  h/a^=l/31,  and  Cf =13 . 685665 

The  situation  with  M=2,  h/ai=l/31,  and  C^=13 . 685665  was 
treated  for  g=0  and  q=c/a^  as  follows: 

1.  For  g=0,  Figure  3  shows  the  comparison  of  the  RAVS 
calculation  (thin  line)  with  Huang's  exact  calculation  (thick 
line)  for  a  step  rise  in  attacking  pressure  with  no  decay  (so 
g=0,  see  (33)).  Note  that  the  curves  are  veiY  close  even  though 
the  minima  of  the  oscillations  in  the  RAVS  calculation  are 
slightly  lower  than  those  in  Huang's  calculation,* 

The  0-th  mode  of  the  radial  displacement  (see  Eq.  5  of 
Huang^)  is  defined  as 


^radial 0  ^ 


2% 

0 


(58) 


because  of  symmetry.  The  nondimens ional  form  of  this  quantity. 


’Because  Huang's  exact  calculation  contains  only  eight  modes, 
it  is  not  quite  exact.  However,  because  the  convergence  is  very 
rapid,  there  is  little  error  in  Huang's  calculation.  Also,  because 
his  curves  were  read  directly  from  his  journal  article,  there  are 
some  small  errors  in  the  presentation  of  his  curves  in  this  paper. 
As  can  be  seen  in  the  program  in  Appendix  C,  the  calculation  with 
the  RAVS  model  had  12  space  intervals  for  6=0  to  0=7C.  For  double 
this  number  of  intervals,  the  calculation  was  about  the  same. 
Thus,  the  calculation  is  sufficiently  accurate  for  this  comparison. 
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_  ^radialO  PC^ 
^radial  0  ~  - - - ~— 

Po 


(59) 


is  plotted  in  Figure  4.  Again,  the  curves  are  very  close.  Also, 
the  RAVS  calculation  at  10  radial  transit  times,  0.1458,  is 
almost  the  same  as  the  exact  asymptotic  value  M/Cf=0 . 1461 
(Huang^)  . 


Figure  3.  Nondimens ional 

velocity  vs.  time, 
0=0,  q=0,  M=2. 


Figure  4.  Nondimens ional 
displacement , 

0-th  mode,  g=D,  M=2 . 


The  1-st  mode  of  the  radial  velocity  (see  Eq.  5  of  Huang*) 
is  defined  as 

1C 

Ui(t)  =  —  J  u(  t,0)  cos0d0  .  (50) 

The  nondimens  ional  form  of  this  quantity, 

Ui  =  ,  (61) 

is  plotted  in  Figure  5.  Again,  the  curves  are  very  close. 

2.  For  g=c/ax,  Figure  6  shows  that  the  agreement  is  also 
very  good  for  q=c/a^. 
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Figure  5.  Nondimensional 

velocity  vs.  time, 
1-st  mode,  g=0,  M=2 . 


Figure  6.  Nondimensional 

velocity  vs.  time, 
6=0,  q=c/a^  ,  M=2 . 


Situation  with  M=4 . 4189 ,  h/ax=l/69 ,  and  €^=12.581197 

The  situation  with  M=4 .4189,  h/a^=l/69 ,  and  C^=12 . 581197  was 
treated  for  g=0  and  g=3c/aj^  as  follows: 

1.  For  g=0,  Figure  7  shows  that  the  agreement  is  very  good 
where  0=0.  In  particular,  the  first 
two  (most  important)  oscillations 
are  reproduced  by  the  RAVS 
calculation  with  very  little 
error.  Also,  the  agreement  at 
Q=iz/2  (Figure  8)  is  very  good.  In 
addition,  the  agreement  at  0=7t  (at 
the  back  of  the  cylinder,  see 
Figure  9)  is  very  good  up  to  two 
radial  transit  times.  For  later  times  the  RAVS  velocity  is  a 
fair  approximation  of  the  exact  velocity  even  though  it  is 
somewhat  lower.  However,  for  most  practical  situations,  the 
front  side  is  of  the  greatest  interest.  For  these  situations, 
the  calculation  for  the  back  of  the  cylinder  does  not  need  to  be 


1.75 

1.5 

1.25 

u 

A 

1 

0.75 

0.5 

/  RAVS 

0.25 

Huang 

T 

2  g— 

S  T(5 

Figure  7.  Nondimensional 

velocity  vs.  time, 
0=0,  g=0,  M=4.4189. 
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Figure  8.  Nondimensional  Figure  9.  Nondimensional 

velocity  vs.  time,  velocity  vs.  time, 

e=7C/2,  q=0,  M=4.4189.  B=n,  q=0,  M=4.4189. 

very  accurate  if  it  is  not  too  large.  In  many  practical 
situations,  the  cylinder  contains  enough  internal  equipment  to 
make  the  entire  object  neutrally  buoyant.  Often,  after  two 
radial  transit  times,  this  internal  equipment  attains  velocities 
that  are  close  to  the  velocity  of  the  shell.  That  is,  the  entire 
object  moves  almost  as  if  it  were  a  neutrally  buoyant  shell  such 
as  was  treated  in  Figures  3-5.  Thus,  if  the  entire  object 
(including  internal  equipment)  were  modeled  using  the  RAVS  model, 
then  it  is  plausible  that  the  motions  would  be  close  to  those  of 
the  neutrally  buoyant  shell  for  which  the  agreement  is  very  good 
with  the  exact  calculations.  Hence,  in  practical  situations,  the 
RAVS  calculation  might  be  close  to  the  exact  motion  even  on  the 
back  of  the  shell . 

The  exact  calculation  of  the  strains  for  this  situation 
(Figures  10-12)  are  in  very  good  agreement  with  the  RAVS 
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calculations  for  all  times  and 
all  locations.  Note  that,  in 
these  figures,  the  ordinate  is 
defined  as  the  strain 
multiplied  by  a  nondimens ional 
factor,  i . e . , 


Figure  10 .  Strain  vs .  non- 
dimensional  time 
e=0,  g=0,  M=4.4189. 


strain* 


strain 


pc^ 

Po 


where 


strain  = 


~  ^radial 


(62) 


Figure  11.  Strain  vs.  non- 

dimensional  time, 
Q=n/2,  q=0,  M=4.4189. 


Figure  12 .  Strain  vs .  non- 

dimensional  time, 
q=0,  M=4.4189. 


The  most  important  oscillatory  behavior  is  accurately  replicated 


*The  small -amplitude  and  high-frequency  oscillations  are  due 
to  using  only  48  intervals  from  0=0  to  0=7i.  If  more  intervals  were 
used,  then  these  amplitudes  would  be  even  smaller  and  these 
frequencies  would  be  even  higher.  This  would  make  the  agreement 
with  the  exact  calculation  even  better. 

Note  that  all  of  the  calculations  have  only  12  intervals 
except  the  strain  calculations  (Figures  10-12  and  Figures  15-17 
which  have  48  intervals)  and  the  calculations  with  q=4.8c/a^ 
(Figures  22-26  which  have  24  intervals) . 
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by  the  RAVS  calculation.  Even  though  the  velocity  (Figure  9) 
at  0=7C  was  somewhat  smaller  than  the  exact  calculation  after  two 
radial  transit  times,  the  strain  (Figure  12)  at  this  location  is 
almost  the  same  as  the  exact  calculation.  This  is  significant 
because  the  strains  are  a  better  measure  of  the  condition  of  the 
shell  than  the  velocity. 

The  0-th  mode  for  the  radial 
displacement  (58)  is  shown  in 
Figure  13 .  Again  the  agreement  is 
very  good.  At  10  radial  transit 
times,  the  RAVS  calculation  is 
0.35111  which  is  almost  the  same 
as  the  exact  asymptotic  value 
(HuangM  ,  M/C^=0 . 35123 .  As 
expected,  the  0-th  mode  of  displacement  is  almost  the  same  as  the 
values  of  strain  for  this  time. 

The  exact  values  of  the  1-th 

mode  for  velocity  (Figure  14)  are 

almost  the  same  as  the  values 

calculated  using  the  RAVS  model  up 

to  two  radial  transit  times.  For 

later  times,  the  RAVS  model  gives 

a  good  approximation  even  though  Figure  14.  Nondimens ional 

velocity  vs.  time, 

the  values  are  somewhat  lower.  1-st  mode, 

q=0,  M=4.4189. 

These  lower  values  are  due  to  the 


Figure  13.  Nondimens ional 
displacement, 
0-th  mode, 
q=0,  M=4.4189. 


lower  values  of  velocity  at  the  back  of  the  cylinder. 
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2.  For  q=3c/ai,  Figures  15-17  show  that  the  agreement  for 
all  locations  and  times  is  also  very  good.  As  in  the  previous 
case,  the  oscillations  are  also  accurately  calculated  by  using 
the  RAVS  method. 


Figure  15.  Strain  vs.  non-  Figure  16.  Strain  vs.  non- 

dimensional  time,  dimensional  time, 

0=0,  q=3c/a^,  M=4.4189.  0=71/2,  q=3c/a^,  M=4.4189. 


Figure  17.  Strain  vs. non- 

dimensional  time, 
0=7t,  q=3c/ai,  M=4 . 4189 . 


Figure  18. 


Nondimens ional 
displacement, 
0-th  mode, 
g=0,  M=6.  41975. 


Situation  with  M=6 .  41975,  h/a^=l/100,  Cf =13 .685665,  and  q=0 

The  situation  with  M=6.  41975,  h/ai=l/100,  Cf =13 . 685665 ,  and 
q=0,  is  compared  for  the  0-th  displacement  mode  in  Figure  18. 
Even  though  the  RAVS  calculation  rises  somewhat  faster  than  the 
exact  calculation,  the  agreement  is  very  good.  The  value  for  10 
radial  transit  times  is  0.46909  which  is  the  same  as  the  exact 
asymptotic  value  (Huang^)  ,  M/Cf =0.46909,  to  five  significant 
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figures.  The  RAVS  calculation  of 
Figure  19,  agrees  very  well  up. to 
At  later  times,  it  is  a  fair 
approximation  even  though  it  is 
lower  than  the  exact  values. 
Unfortunately,  Huang*  did  not 
present  the  calculations  for  the 
velocities  and  strains  at 
individual  locations  for  this 
situation.  As  previously 


the  1-st  velocity  mode, 
two  radial  transit  times. 


Figure  19.  Nondimens ional 

velocity  vs .  time , 
1-st  mode, 
q=0,  M=6. 41975. 


explained,  for  many  situations  of 

interest  with  internal  equipment  inside  the  shell,  the  velocities 
at  later  times  would  tend  to  be  the  same  as  in  the  neutrally 
buoyant  situation.  For  the  neutrally  buoyant  situation,  the 
RAVS  calculation  agrees  very  well  with  the  exact  calculation. 
Thus,  the  small  deviation  after  two  radial  transit  times  might  be 
of  little  practical  interest.  Also,  in  the  previous  situation, 
the  velocities  on  the  front  of  the  cylinder  and  the  strains  at 
all  locations  were  very  close  to  the  exact  calculations,  even 
though  the  0-th  velocity  mode  was  somewhat  lower  than  the  exact 
calculation  after  two  radial  transit  times.  Thus,  it  is 
plausible  that  the  velocities  on  the  front  of  the  cylinder  and 
the  strains  at  all  locations  are  in  better  agreement  than  the  1- 
st  velocity  mode. 


Situation  with  M=9 . 09184 ,  h/ai=l/142 ,  €^=12.581197,  and  g=0 
Similar  remarks  apply  to  the  situation  with  M=9 .09184, 
h/ai=l/142,  Cf =12 . 581197,  and  q=0,  which  is  presented  in  Figures 
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20  and  21.  In  Figure  20,  the  value  for  10  radial  transit  times 
is  0.7213  which  is  almost  the  same  as  the  exact  asymptotic  value 
(Huang'*)  ,  M/(f=0.  7227. 


Figure  20.  Nondimens ional 
displacement , 
0-th  mode, 
q=0,  M=9. 09184. 


Figure  21.  Nondimens ional 

velocity  vs.  time 
1-st  mode, 
q=0,  M=9. 09184. 


Situation  with  M=4 . 4189 ,  h/ax=l/69 ,  C^=12 .581197 ,  and  q=4 . 8c/a^ 

The  velocities  for  the  situation  with  M=4 . 4189,  h/ax=l/69 , 
€^=12.581197,  and  q=4 . 8c/a^  were  calculated  for  various  locations 
and  presented  in  Figures  22-26.  These  figures  show  that  the  RAVS 
calculations  agree  very  well  with  the  exact  calculations  for  all 
locations  and  times. 


Figure  22.  Nondimens ional 

velocity  vs.  time, 
0=0,  q=4.8c/a^,  M=4.4189. 


Figure  23.  Nondimens ional 

velocity  vs.  time, 
6=7t/6,  q=4.8c/a^,  M=4.4189. 
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Figure  24.  Nondimens ional 

velocity  vs.  time, 
0=7t/3,  q=4.8c/a^,  M=4.4189. 


Figure  25.  Nondimens ional 

velocity  vs.  time, 
0=271/3,  q=4.8c/a,,  M=4.4189. 


Note  that  the  RAVS  calculation  for 
this  situation  had  24  space 
intervals  from  0=0  to  0=7i  for  the 
finite-difference  calculation  as 
programmed  in  Appendix  C. 


Figure  26.  Nondimens ional 

velocity  vs.  time, 

0=7t,  q=4.8c/a^,  M=4.4189. 

RIGID  CYLINDERS 

It  is  shown  in  Ref.  1  that  the  RAVS  calculations  for  rigid 
(and  immovable)  cylinders  attacked  by  plane  waves  are  in  good 
agreement  with  the  exact  calculations.  However,  the  preceding 
calculations  for  elastic  cylinders  attacked  by  planar  waves  are 
in  even  better  agreement  with  the  exact  calculations.  Also,  in 
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Ref.  1,  RAVS  calculations  (dashed  lines)  are  compared  to 
Huang's"'  exact  calculations  (solid  lines)  for  a  spherical  wave 
from  a  source  that  is  0.5  from  the  nearest  point  on  the  cylinder 
(see  (2)  with  Phydro=Prad=0 ,  a^=l,  and  c=I)  and  has  a  step  pressure 
p(t)=0.5/r,  for  Ost,  see  Figure  27.  It  is  seen  that  the  error 
for  0  degrees  is  very  small.  The  exact  pressures  at  90,  120,  and 
180  degrees  are  so  much  smaller  than  that  at  0  degrees,  that  for 
many  purposes  they  can  be  disregarded  (at  these  angles,  the  RAVS 
calculations  give  zero  pressure) .  Similar  remarks  apply  to 
Figure  28.  Because  the  RAVS  calculations  for  the  elastic  case 
are  in  better  agreement  with  the  exact  calculations  for  planar 
waves,  it  is  plausible  that  the  agreement  for  the  elastic  case 
would  be  in  better  agreement  for  spherical  waves.  Thus,  it  is 
plausible  that  the  RAVS  calculations  for  spherical  waves  would  be 
almost  the  same  as  the  exact  calculations.  This  is  important 
because  spherical  waves  are  of  much  more  practical  interest  than 
planar  waves. 


Figure  27 .  Rigid  cylinder 

attacked  by  a  spherical  wave 
with  constant  pressure, 
nearest  distance  =  0.5aj, 
longitudinal  distance  =  0 . 


Figure  28.  Rigid  cylinder 

attacked  by  a  spherical  wave 
with  constant  pressure, 
nearest  distance  =  O.Saj, 
longitudinal  distance  =  1.41421aj. 
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SUMMARY  AND  CONCLUSIONS 


An  equation  of  motion  (1)  for  a  structural  surface  was 
presented.  The  equations  for  the  pressure  that  develops  when  an 
acoustic  wave  interacts  with  a  curved  and  compliant  surface  were 
presented  ((2) -(14)).  These  equations  were  derived  using  the 
RAVS  model. With  these  equations,  an  expression  for  the 
velocity  ((15) -(20))  of  the  structural  surface  was  derived  from 
the  equation  of  motion,  and  a  numerical  method  was  presented  for 
using  this  expression.  This  was  applied  to  the  case  of  a  planar 
wave  attacking  an  elastic  cylindrical  shell.  Hence,  a  computer 
program  was  written  for  this  problem.  Calculations  from  this 
program  (using  the  RAVS  model)  were  shown  to  be  in  good  agreement 
with  all  of  the  exact  calculations  presented  by  Huang*  for  all 
locations  and  times.  Thus,  the  RAVS  model  gives  a  good 
approximation  for  an  elastic  cylinder  being  attacked  by  a  planar 
wave . 

It  was  shown  that  the  RAVS  approximation  is  even  better  for 
spherical  attacking  waves  which  are  of  much  more  practical 
interest.  In  these  situations,  the  source  of  the  attacking  wave 
is  a  finite  distance  from  the  structural  surface.  In  the  limit 
of  the  source  being  very  close  to  a  structural  surface,  the 
structural  surface  can  be  considered  as  planar.  In  this 
situation,  the  RAVS  approximation  becomes  exact  for  both  a  rigid 
surface  and  a  free  surface. 

Of  course,  the  results  of  the  RAVS  model  should  be  compared 
to  other  exact  calculations  such  as  those  for  an  elastic 
spherical  shell . ® 
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APPENDIX  A 


SOLUTION  OF  THE  EQUATION  OF  MOTION 


Substituting  (14)  (for  Pj-ad)  in  (2)  (for  Ptot)  and  the 
resulting  expression  in  the  equation  of  motion  (1)  gives 


m 


uit)  =Pt,ydzo^pU)  -e  ^  /  e  *  ^cosp^^  dt>-p^^^{t)  . 


-r/c 


(63) 


Although,  cosp  in  (63)  must  be  positive,  it  might  be  small  and 
cause  difficulties  (such  as  numerical  instabilities)  in  numerical 
calculations  because  it  appears  with  the  unknown  velocity  on  the 
right-hand  side  of  this  equation.  To  avoid  these  difficulties, 
the  equation  of  motion  for  the  structural  surface  is  expressed  in 
an  equivalent  form  that  is  well-behaved  even  when  cosp  is  small. 
To  do  this,  (63)  is  transformed  so  that  u(t)  does  not  appear  on 
the  right-hand  side  of  the  equation.  The  contribution  of  the 
pressure  of  the  scattered  wave  in  this  equation  is  then  expressed 
as  a  function  of  the  attacking  pressure. 


TRANSFORMED  EQUATION  OF  MOTION 

Equation  (63)  can  be  written  as 


Ct 

e  *  u(t) 


Pg(fc) 

m 


pcujt') 

mcosp 


dt'  , 


(64) 


where 


Pgit)  =  pit)  +p^^at(t)  -  [Pscr<t)  -Phydro'i  .  (65) 

Note  that  the  bracketed  term  in  this  equation  is  zero  until  the 
attacking  wave  arrives  at  the  structure.  Differentiating  both 
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sides  of  (64)  gives 


_d 

dt 


Ct 

^  u{t) 


incosP 


_d 

dt 


P,(t) 


m 


(66) 


Solving  this  equation  gives 


Ct 


-  -P  g  - 


e  *  u(t)  =  e  f  e 

-rtc 


d 

mcosp  ^ 

dV 


m 


dt* 


m  m  cosP 


L 


m 


(67) 


Thus,  we  obtain  a  transformed  equation  of  motion, 


c 

mu(t)  =  Pg(t)  -  (t)  -  j*  e^^'pg(t')  dt' 

'  '  -r/c 


(68) 


where  (see  (65)) 

t 

J^(t)  =©“’>"/  e’>^p^(tOdt' =  J(t)  +J^^,,(t)  -I^,^(t)  ,  (69) 

-r/e 


in  which 


t 

J(t)  =  e^^'pit')  dt'  ,  (70) 

0 


t 

I,,,,{t)  =  e’'*^  p^^,,(t')  dt'  ,  (71) 

0 


t 

-r/c 


and 
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(73) 


n 

h  incosp 

is  called  "the  damping  coefficient."*  The  lower  limit  of  the 
integration  in  (70)  and  (71)  is  zero  because  the  pressure  of  the 
attacking  wave  and,  thus,  the  pressure  of  the  scattered  wave  are 
zero  before  their  arrival  at  t=0.  Therefore,  the  transformed 
equation  of  motion  (68)  can  be  written  as 

milit)  =j(t)+|j(t)+4^,e(t)+|^scat(t)-4tr(t)-|j^tr(t)  .  (74) 

SCATTERED  PRESSURE  CONTRIBUTION  AS  A  FUNCTION 
OF  ATTACKING  PRESSURE 

From  (4) ,  it  is  seen  that 

.  (75) 

By  performing  the  integration, 

)dt',  (76) 

0 

on  all  of  the  terms  in  (75) ,  it  is  seen  that 

*It  can  be  seen  from  (69)  that 

4(t)  +TiJ^(t)  =  p^{t)  , 

which  shows  that  T|  behaves  as  a  damping  coefficient. 
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*The  pressure  of  the  incident  wave  and,  thus,  the  pressure 
of  the  scattered  wave  are  zero  in  the  limit  of  positive  time  going 
to  zero.  If  there  is  a  sudden  rise  in  the  pressure  at  time  zero, 
the  pressure  can  be  represented  (to  any  desired  accuracy)  by  a 
function  that  is  zero  in  the  limit  of  positive  time  going  to  zero. 
The  function  would  then  increase  at  an  arbitrarily  high  rate  of 
change  up  to  its  first  peak  pressure. 
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EQUIVALENT  EQUATION  OF  MOTION 

Substituting  (80)  in  (74)  gives  an  equivalent  equation  of 
motion, 

mu(t)  =  j(t) +-2j{t) +[r(t)+-J(t)l-2^-j^t^(t) ,  (81) 
b  [  r  J  cosp  b 

where  I(t)  and  determined  by  (70)  and  (72).  The 

right-hand  side  of  (81)  does  not  involve  the  (unknown)  velocity 
(as  in  (63)).  However,  the  right-hand  side  does  require  the  mass 
density  of  the  structural  surface  and  the  structural  resistance 
pressure.  Integration  of  both  sides  of  (81)  and  dividing  both 
sides  of  the  resulting  expression  by  m  gives  the  velocity  of  the 
structural  surface. 


ii(t)  = 


Kt)  /i. 

cos6 

m  \ 

cosp 

m 

mb 

-r/c  * 


see  (15) -(20)  of  the  main  text. 
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APPENDIX  B 


STRUCTURAL  RESISTANCE  PRESSURE 


In  this  appendix,  the  notation  in  Huang's  paper*  is  used 
for  the  expressions  in  the  first  section.  In  the  second  section, 
these  expressions  are  converted  to  the  notation  of  the  main  text 
to  obtain  the  structural  resistance  pressure, 

HUANG'S  NOTATION 

Equations  (8)  and  (9)  of  Huang's  paper*  state  that 


dT^ 


dT^ 


Cn^n 


D  n 


=  0 


n  =  0,1,2, .. . 
where 

A„  =  cMl+u^hVOa*)] 

Bjj  =  Djj  =  nC^  [1  +n^h^/  Oa^)  ] 

C„  =  +iiV(3a2)]  , 


(83) 


where  Eq.  (5)  of  Huang's  paper  states  that 


P(l,0,r)  =  P„(1,D  cos(ne)  , 

i2=0 

•• 

w(0,D  =  ]^  w„(D  cos(j30)  , 

a»0 

«• 

v(0,r)  =  vi,(r)  sin(n0)  , 
and  Eq.  (1)  of  Huang's  paper  states  that 


(84) 


33 


T  =  ct/a  ,  w  =  hja 


v=hja,  P(l,e,r)  =  p(a,0,  t)  /  (pc^)  , 


M=pa/(2hp^),  C2  =  p/ [p^(l -v2)  c2]  ,  (85) 


where 


c  =  acoustic  wave  speed, 

t  =  time  after  arrival  of  the  wave  front  at  the  shell, 
a  =  middle  surface  radius  of  the  cylindrical  shell, 

=  radial  deflection  , 

8g  =  circumferential  deflection, 
p(a,0,  t)  =  total  pressure  in  the  fluid  at  the  shell, 
p  =  unperturbed  mass  density  of  the  fluid, 
h  =  half  thickness  of  the  shell, 

Pg  =  mass  density  of  the  shell  material, 

E  =  Young's  modulus  of  the  shell  material,  and 
V  =  Poisson's  ratio  of  shell  material. 

(86) 


Performing  the  operation, 


^  cos (n0)  ,  (87) 

n-0 

on  both  sides  of  the  first  equation  in  (83)  and  using  (84)  yields 

+  c^w  +  =  -AfP(l,0,r)  ,  (88) 

dT^  3a^  3a^ 

where  the  superscripts  (in  parentheses)  indicate  partial 

differentiation  with  respect  to  6.  Also,  performing  the 

operation. 
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(89) 


5^  sin(j20)  , 

n-1 

on  both  sides  of  the  second  equation  in  (83)  yields 

dT^  3a2j  3a2 

Substituting  (85),  in  (88)  and  (90)  and  multiplying  the 
equations  by 

a/c^ 

gives 

fir  +  -^fir‘^^  +  fie'^^  -  — fie‘^^1  = 

.  ^  3a2  ®  3a2  ®  J 

and 


2hp, 


at2 


2Eh 


(1  -v2)a 


^fie 

2hp^— I 
®  at2 


2Eh 


(1  -  v^)  a^ 


H-#) 


3a 


.5^(3) 


CONVERSION  TO  NOTATION  IN  MAIN  TEXT 

The  following  variables  in  the  main  text  [see  (1)  , 
Ptoti-  see  (48),  for  and  dg,-  and  see  (7)  for  aj  ,  are 

defined  using  Huang's  notation: 

m  =  2hpg  ,  Huang, 

^radial  =  “fir  ^  Huang, 

=  Sg  ,  Huang 
a^  =  a  ,  Huang, 

Ptot  =  P  /  Huang. 


(90) 

resulting 

(91) 

-p  (92) 

0  .  (93) 

for  m  and 
now 

(94) 
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Thus,  (92)  becomes 


-m- 


d^d 


radial 


2  Eh 


dt‘ 


(1  -  v^) 


-d 


radial 


j(4) 

^radial 


3ai 


(3) 


Ptot 

(95) 


or 


_  d^d^adial  _  _  _ 

^  Q^2  ~  ~  Pstr  > 


(96) 


where  (see  the  equation  of  motion  (1)) 


2Eh 


'stT 


(l-v2)  aj^ 


^radial  2 

3  a? 


^radial 


3  a^^ 


(3) 


(97) 


is  the  structural  resistance  pressure,  see  (48)  in  the  main  text, 
Also,  (93)  becomes 


d^cL 

m — —  + 


2  Eh 


dt^  (l-v^)ai^ 


lrf<2)  + 

pe  +  ^radial  — - X  ^radial 

I  3ai2 


=  0  ,  (98) 


see  (49)  in  the  main  text. 
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APPENDIX  C 


COMPUTER  PROGRAM 


The  following  finite-difference  program  is  written  in  the 
Mathematical  language.  Because  the  variable  names  are  similar  to 
those  in  the  text,  it  should  be  fairly  easy  to  follow  and  modify. 
A  sample  run  is  also  included. 


Date [ ] 

Print [ " " ] ; 

Print  [  '■ 

THE^  VIRTUAL-SOURCE 

MODEL 

"]  ; 

Print  [  '■  "  ]  ; 

Print [ " 

APPLIED  TO  AN  ELASTIC 

CYLINDER 

"]  ; 

Print [ " " ] ; 

Print [ ” 

FOR  A  PLANAR  SHOCK 

WAVE 

"] ; 

Print [ " " ] ; 

(*  See  Huang,  Eq.  1  and  his  discussion  before  Eq.  25.  *) 
(*  MHuang=rho*al/ (2*h*rhos) ; 

CHuang2=E/  (rhos*  (l-nu''2)  *c''2) )  *) 

(*  MHuang=2;  ha=l/31;  CHuang2=13 . 685665;  *) 

MHuang=4 . 41890 ;  ha=l/69;  CHuang2=12 . 581197 ; 

(*  MHuang=6 . 41975 ;  ha=l/100;  CHuang2=13 . 685665; *) 

(*  MHuang=9. 09184;  ha=l/142;  CHuang2=12 . 581197 ;* ) 
al=l.;  rho=l.;  c=l.;  (*  Non-dimensional  results  are  the 

same  for  all  positive  values  of 
al,  rho,  and  c.  *) 
q=0;  (*  q=c/al;  3*c/al;  4.8*c/al; 

q  is  the  same  as  beta*c/al  in  Huang,  Eq.  25.  *) 

Omega=l.;  (*  Non-dimensional  results  are  the  same  for  all 
non- zero  values  of  Omega.  *) 
pO=Omega;  (*  See  Huang,  Eq.  25.  *) 

mass=rho*al/MHuang;  kH= {rho*c*c/al) *CHuang2 /MHuang; 

rhocmass=rho*c/mass; 

h2=ha*ha/3 ; 

nmax=12;  (*  =  maximum  number  of  space  intervals 
from  0  to  Pi  *) 

Dtheta=N[Pi] /nmax; 

mmax=750;  (*  =  maximum  number  of  time  intervals 

from  0  to  ttmax  *) 

ttmax=:2 . 5*al/c; 

D t ime= t tmax /mmax ; 

mPrint=15;  (*  =  interval  between  printing  *) 
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DyDx  [y_,  Deltax_,  sgn_]  :  =Module  [  {dummy}  , 

(*  sgn=l,  y  even;  sgn=-l,  y  odd  *) 
ya=Drop  [  Join  [sgn*Take  [y,  {2 , 2}  ]  ,y]  ,  -1]  ; 
Lml=Length [y] -1 ; 

yb=Drop  [Join  [y,  sgn*Take  [y,  (Lml,  Lml}  ]  ]  ,  1]  ; 
(yb-ya)  /  (2*Deltax)  (*  =  rate  of  change  of  y  *)] 

Sh[x_]  :=If [x==0,  1,  Sinh[x]/x  ] 

etaF  [  theta_]  :  =Module  [  (dummy) , 

S2= (Sin [theta/2] ) ^2;  SS=Sqrt [1+8 . *S2 ] ; 
b=al*SS;  cbeta= (1+2*S2) /SS; 

N[  (c/b)  +  (rho*c/ (mass*cbeta) ) ] ] 


cbF  [  theta_]  :  =Module  [  {dummy} , 

S2=( Sin [theta/2] )^2;  SS=Sqrt [1+8 . *S2 ] ; 
b=al*SS; 

N[c/b]] 

edF  [  theta_]  :  ^Module  [  {dummy} , 
eta=etaF [theta] ; 

2*Exp  [-eta*Dtime/2]  *Sinh[eta*Dtime/2]  /eta] 

JOF  [  t_]  ;  =Module  [  {dummy}  , 

If[chi''2  >  (epsilon*s)  ^2 , 
eq=Exp [-q*t/2] *Sh[q*t/2] ; 
ee=Exp[-eta*t/2]Sh[eta*t/2] ; 

J0=p0* (eq-ee) *t/ (eta-q) ,  J0=p0*t*t/2]  ] 

uOF[theta_,  tt_]  :  =Module  [  {dummy} , 

S2=  (Sin [theta/2 ]) ''2;  rMrO=2*al*S2 ;  t=tt-rMrO/c; 

I  f  [  t  <  0 ,  0 , 

SS=Sqrt [1+8*S2]  ;  b=al*SS;  cbeta= (1+2*S2) /SS; 
eta=  (c/b)  +  (rho*c/  (mass*cbeta) )  ; 
s=(eta+q)/2;  chi= (eta-q) /2  ; 

IO=pO*t*Exp[-s*t] *Sh[chi*t]  ; 
epsilon=0 . 0001 ; 

If[q*t  >  epsilon, 

If[eta*t  >  epsilon, 

cs=Cosh[chi*t] +s*t*Sh[chi*t] ; 

J0=p0* (l-Exp[-s*t] *cs) / (q*eta) , J0F[t] ]  ,  J0F[t]  ] 
ctheta=Cos [theta]  ; 
u0m=l0*  ( 1+ctheta/cbeta)  +  (c/b)  *J0; 
uOm/mass  ]  ] 

etaT=Table  [etaF  [n*Dtheta]  ,  {n,0,nmax}]  ; 
edT=Table  [edF  [n*Dtheta]  ,  {n,0,nmax}]  ; 
cbT=Table  [cbF  [n*Dtheta]  ,  {n,0,nraax}]  ; 

u=Table[0,  {n,  0,nmax}]  ;  ustr=u;  Istr=u;  dr=u;  v=u;  dth=u 
urc0=strainrc20= 

Table  [0  .  ,  {n,  0 , nmax} ,  {m,  0, Floor  [mmax/mPrint] }  ]  ; 
dr 0rc2Oa=ulrcO=Table  [  0  . ,  {m,  0 ,  Floor  [mmax/mPr int  ]}]; 

Do [tt=N[m*Dtime] ; 
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uO=Table [uOF [n*Dtheta, tt] , {n,0,nmax}]  ; 
drl=DyDx [dr , Dtheta,  1]  ;  dr2=DyDx[drl, Dtheta, -1]  ; 
dr3=DyDx[dr2,Dtheta,  1]  ;  dr4=DyDx  [dr3  ,  Dtheta, -1]  ; 
dthl=DyDx [dth, Dtheta, -1]  ;  dth2=DyDx [dthl, Dtheta,  1]  ; 
dth3=DyDx[dth2,Dtheta,-li ; 

pstr=kH* (dr+h2*dr4-dthl+h2*dth3) ;  strain= (dthl-dr) /al; 
DIstr= (pstr-etaT*Istr) *edT;  Istr=Istr+DIstr ; 

Dustr= (DIstr+Istr*cbT*Dtime) /mass; 
ustr=ustr+Dustr ;  uL=u;  u=uO-ustr; 

If [ IntegerQ [m/mPrint ] ,  mP=m/mPrint ; 

Print [ "  "];  Print [ "c*tt/al  =  ",c*tt/al]; 

Print  [  "u*rho*c/Omega  =  ",Chop[u*rho*c/C)mega]]; 
Print  [  "strain*rho*c''2 /Omega  =  ", 

Chop  [strain*rho*c''2 /Omega]  ]  ; 
drO=( .5*(dr[ [1] ]+dr[ [nmax+1] ] ) 

+Siim[dr  [  [n]  ]  ,  {n,  2,nmax}  ] )  /nmax; 
ul=( .5* {u[ [1] ] -u[ [nmax+1] ] ) 

+Siim[u[  [n]  ]  *Cos  [  (n-l)  *N[Pi]  /nmax] 

, {n, 2 , nmax} ] ) *2 /nmax; 

Print [ 

"dr0*rho*c^2/ (Omega*al)  =  ", 
dr0*rho*c^2/ (Omega*al)  , 

"  ul*rho*c/Omega=  “  ,ul*rho*c/Omega] ; 

Do[  urcO[ [n+l,mP+l] ] 

=u [ [n+1] ] *rho*c/Omega; 
strainrc20 [ [n+l,mP+l] ] 

=s train  [  [n+l]  ]  *rho*c'^2 /Omega; 

, {n, 0 , nmax) ] ; 

drOrc20a  [  [mP+1]  ]  =dr0*rho*c''2/  (Omega*al)  ; 
ulrcO[ [mP+1] ] =ui*rho*c /Omega; 

, dumb=0] ; 

dr=dr+ (u+uL) *Dtime/2 ; 

vd=(kH/mass) * ( (l+h2) *dth2-drl+h2*dr3)  ; 

vL=v;  v=v+vd*Dtime;  dth=dth+ (v+vL)  *Dtime/2  ; 

,  {m, 0,mmax} ]  ; 
urcO  »  urc0.dat; 
strainrc20  »  strc20.dat; 
drOrc20a  »  drc20a.dat; 
ulrcO  »  ulrc0.dat; 

Print ["DONE"] 


{1996,  1,  30,  16,  59,  28} 

THE  VIRTUAL-SOURCE  MODEL 
APPLIED  TO  AN  ELASTIC  CYLINDER 
FOR  A  PLANAR  SHOCK  WAVE 


C*tt/al  =  0 

u*rho*c/Qmega  =  {0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0} 
strain*rho*c''2 /Omega  =  {0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0} 
dr0*rho*c''2/ (Omega*al)  =  0  ul*rho*c/Omega=  0 
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c*tt/al 


0.05 


-6  -7 

u*rho*c/Omega  =  {0.390278,  0.134938,  -1.754  10  ,  1.31981  10  , 

-6  -8 

-1.73709  10  ,  -4.87924  10  ,  0,  0,  0,  0,  0,  0,  0} 

strain*rho*c^2 /Omega  =  {-0.00878012,  -0.000692333,  -0.0000658705, 

-7  -7  -10  -10 

-2.6909  10  ,  -1.34421  10  ,  1.11087  10  ,  -1.11216  10  ,  0,  0,  0, 

0,  0,  0} 

dr0*rho*c''2/ (Omega*al)  =  0.000429056  ul*rho*c/Omega=  0.0542462 
c*tt/al  =  0.1 

u*rho*c/Omega  =  {0.68735,  0.491574,  -0.00024465,  -0.0000248791, 

-6  -8  -10 

-0.000015459,  -4.00469  10  ,  -1.20485  10  ,  -1.84257  10  ,  0,  0,  0, 

0,  0} 

strain*rho*c^2/Omega  =  {-0.0327916,  -0.0153423,  -0.00119008, 

-7  -8 

-0.000216349,  -0.0000148417,  -9.67291  10  ,  -9.06537  10  , 

-9  -10 

-2.01236  10  ,  -3.1317  10  ,  0,  0,  0,  0} 

dr0*rho*c^2/ (Omega*al)  =  0.00276337  ul*rho*c/Omega=  0.136377 

c*tt/al  =  0.15 

u*rho*c/Cttnega  =  {0.906989,  0.758442,  0.130266,  -0.000555736,  -0.00008591, 

-7  -8  -9 

-0.0000269907,  -6.22903  10  ,  -5.07746  10  ,  -3.78493  10  ,  0,  0,  0,  0} 

strain*rho*c^2 /Omega  =  {-0.0627982,  -0.0440431,  -0.00636443,  -0.00216181, 

-6  -7 

-0.000180932,  -0.0000378213,  -2.85768  10  ,  -3.28619  10  , 

-8  -9  -10 

-2.63074  10  ,  -1.64085  10  ,  -1.55108  10  ,  0,  0} 

dr0*rho*c''2/ (Omega*al)  =  0.00701587  ul*rho*c/Omega=  0.216411 

c*tt/al  =  0.2 

u*rho*c/Omega  =  {1.06626,  0.950778,  0.469704,  -0.00322633,  -0.000435931, 

-6  -7  -8 

-0.000146382,  -0.000011791,  -1.47  10  ,  -1.13404  10  ,  -1.31009  10  , 

-9 

-1.1139  10  ,  0,  0} 

strain*rho*c''2 /Omega  =  {-0.0900519,  -0.0796894,  -0.0304906,  -0.00849706, 

-6 

-0.00117721,  -0.000334082,  -0.0000303624,  -6.61006  10  , 

-7  -8  -9  -10 

-5.26616  10  ,  -7.67038  10  ,  -6.06395  10  ,  -5.77085  10  ,  0} 

dr0*rho*c''2/ {Omega*al)  =  0.013771  ul*rho*c/Omega=  0.309292 

c*tt/al  =  0.25 

u*rho*c/Omega  =  {1.18173,  1.08383,  0.713623,  -0.0108063,  -0.00207721, 

-6 

-0.000665031,  -0.0000820345,  -0.0000158027,  -1.40114  10  , 

-7  -8  -9  -10 

-2.7135  10  ,  -2.21772  10  ,  -3.20792  10  ,  -5.12757  10  } 

strain*rho*c''2 /Omega  =  {-0.110981,  -0.115601,  -0.0727123,  -0.021653, 

-6 

-0.00536965,  -0.00153255,  -0.00020353,  -0.0000547712,  -5.25935  10  , 

-6  -8  -8  -9 

-1.15562  10  ,  -9.65689  10  ,  -1.60209  10  ,  -2.54973  10  } 

dr0*rho*c''2/ {Omega*al)  =  0.0227186  ul*rho*c/Omega=  0.374488 

c*tt/al  =0.3 

u*rho*c/Omega  =  {1.26616,  1.17245,  0.877804,  0.029222,  -0.00758226, 

-6 

-0.00234143,  -0.000401832,  -0.0000971238,  -0.0000110074,  -2.6893  10  , 

-7  -8  -9 

-2.49546  10  ,  -5.09881  10  ,  -8.47793  10  } 

strain*rho*c^2/Ctaiega  =  {-0.128169,  -0.146546,  -0.122157,  -0.0426228, 
-0.016362,  -0.00480748,  -0.00096836,  -0.000274106,  -0.0000354739, 

-6  -7  -7  -8 
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-9.23502  10  ,  -9.20385  10  ,  -2.06215  10  ,  -3.53533  10  } 

dr0*rho*c''2/ (Omega*al)  =  0.0331557  ul*rho*c/Omega=  0.42368 

c*tt/al  =  0.35 

u*rho*c /Omega  =  {1.32684,  1.23031,  0.980321,  0.338326,  -0.020849, 

-0.00655868,  -0.00154461,  -0.000417103,  -0.000062901,  -0.0000169489, 

-6  -7  -8 

-1.93476  10  ,  -4.75978  10  ,  -8.93905  10  } 

strain*rho*c^2 /Omega  =  {-0.147291,  -0.169524,  -0.167398,  -0.0802541, 
-0.0370165,  -0.0118034,  -0.00337508,  -0.000985651,  -0.000176098, 

-6  -6  -7 

-0.0000489309,  -6.22838  10  ,  -1.62146  10  ,  -3.24261  10  } 

dr0*rho*c'"2/ (Omega*al)  =  0.0453529  ul*rho*c/Omega=  0.48801 

c*tt/al  =0.4 

u*rho*c/Omega  =  {1.36622,  1.26871,  1.04042,  0.549674,  -0.0458149, 
-0.0155475,  -0.00475985,  -0.001385,  -0.000275576,  -0.0000779002, 

-6  -7 

-0.0000112943,  -3.06547  10  ,  -6.83805  10  } 

strain*rho*c''2 /Omega  =  {-0.173062,  -0.184802,  -0.200042,  -0.133486, 
-0.0676164,  -0.0250382,  -0.00914231,  -0.00281387,  -0.000671982, 

-6  -6 

-0.000193887,  -0.0000320787,  -9.02277  10  ,  -2.1991  10  } 

dr0*rho*c^2/ (Omega*al)  =  0.0591983  ul*rho*c/Omega=  0.528656 

c*tt/al  =  0.45 

u*rho*c/Omega  =  {1.38455,  1.29527,  1.07588,  0.681945,  -0.0847244, 
-0.0326425,-0.0121164,  -0.00380232,  -0.000964665,  -0.000283381, 

-6 

-0.0000520526,  -0.0000149929,  -4.05749  10  } 

strain*rho*c"2 /Omega  =  {-0.206103,  -0.195723,  -0.2172,  -0.189263, 
-0.105065,  -0.0473296,  -0.0203005,  -0.00681753,  -0.00205813, 
-0.000619582,  -0.000131101,  -0.0000388247,  -0.0000116707} 
dr0*rho*c^2/ (Omega*al)  =  0.0739668  ul*rho*c/Omega=  0.551379 

c*tt/al  =  0.5 

u*rho*c/Omega  =  {1.38303,  1.31357,  1.10023,  0.756233,  -0.136777, 
-0.0616435,  -0.0262858,  -0.00905239,  -0.00279995,  -0.000861915, 
-0.000196245,  -0.0000592985,  -0.0000194509} 
strain*rho*c^2 /Omega  =  {-0.24203,  -0.206774,  -0.221796,  -0.23533, 
-0.143549,  -0.0798836,  -0.0384272,  -0.0145521,  -0.00525634, 
-0.00168049,  -0.000440301,  -0.000137011,  -0.000050294} 
dr0*rho*c^2/ (Omega*al)  =  0.0890722  ul*rho*c/Omega=  0.561356 

c*tt/al  =  0.55 

u*rho*c/Omega  =  {1.36612,  1.32405,  1.121,  0.794122,  0.107007,  -0.105335, 
-0.0498887,  -0.019213,  -0.00694996,  -0.00227573,  -0.00062361, 
-0.000198169,  -0.0000776547} 

strain*rho*c''2/Omega  =  {-0.273095,  -0.221529,  -0.22066,  -0.263712, 

-0.183344,  -0.120828,  -0.063793,  -0.0278615,  -0.0115378,  -0.00400039, 
-0.00125217,  -0.000413182,  -0.000181248} 
dr0*rho*c''2/ {Omega*al)  =  0.104631  ul*rho*c/Omega=  0.588568 

c*tt/al  =  0.6 

u*rho*c/Omega  =  {1.34179,  1.32567,  1.13979,  0.814657,  0.270561, 

-0.163875,  -0.084728,  -0.0369165,  -0.0151172,  -0.00534854,  -0.0017124, 
-0.00057744,  -0.000264697} 

s train* rho*c'' 2 /Omega  =  {-0.291559,  -0.241203,  -0.221837,  -0.272703, 

-0.22432,  -0.164998,  -0.0954124,  -0.0483004,  -0.0222821,  -0.00853444, 
-0.00309035,  -0.001096,  -0.00055979} 
dr0*rho*c''2/ (Omega*al)  =  0.12082  ul*rho*c/Omega=  0.605078 

c*tt/al  =  0.65 

u*rho*c/Omega  =  {1.31976,  1.31769,  1.1537,  0.831515,  0.375294,  -0.233867, 
-0.131439,  -0.0648301,  -0.0293742,  -0.0113763,  -0.00414687, 
-0.00150054,  -0.00078646} 

strain*rho*c-'2 /Omega  =  {-0.293339,  -0.264247,  -0.23159,  -0.266693, 

-0.257485,  -0.205179,  -0.131377,  -0.0763083,  -0.0385975,  -0.0165318, 
-0.00675127,  -0.00260922,  -0.00151212} 
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dr0*rho*c^2/ (Omega*al)  =  0.137066  ul*rho*c/Omega=  0.615307 


c*tt/al  =  0.7 

u*rho*c/Omega  =  {1.30839,  1.30102,  1.15827,  0.851501,  0.441232, 

-0.308667,  -0.189149,  -0.104791,  -0.0518179,  -0.0221458,  -0.00900568, 
-0.00353419,  -0.00207205} 

strain*rho*c''2 /Omega  =  (-0.280369,  -0.287027,  -0.251465,  -0.253976, 
-0.277928,  -0.234403,  -0.168588,  -0.110467,  -0.060931,  -0.0293271, 
-0.0132694,  -0.00565358,  -0.00363233} 
dr0*rho*c''2/ (Omega*al)  =  0.15292  ul*rho*c/Omega=  0.622956 

c*tt/al  =  0.75 

u*rho*c/Omega  =  (1.31148,  1.27865,  1.15075,  0.874968,  0.484809, 

-0.330613,  -0.25518,  -0.15684,  -0.0840954,  -0.0397728,  -0.0177835, 
-0.00763362,  -0.00491012} 

strain*rho*c''2 /Omega  =  (-0.260442,  -0.305222,  -0.277531,  -0.243406, 
-0.285732,  -0.248481,  -0.203013,  -0.147248,  -0.0887421,  -0.047942, 
-0.0237822,  -0.0112609,  -0.0078662} 
dr0*rho*c''2/ (Omega*al)  =  0.168043  ul*rho*c/Omega=  0.633154 

c*tt/al  =  0.8 

u*rho*c/Omega  =  (1.32667,  1.25515,  1.1324,  0.897965,  0.517252,  -0.151369, 
-0.325165,  -0.218549,  -0.126886,  -0.0663104,  -0.0323033,  -0.0152529, 
-0.0105897} 

strain*rho*c^2 /Omega  =  (-0.244415,  -0.315389,  -0.302186,  -0.241256, 
-0.284926,  -0.255,  -0.230604,  -0.181612,  -0.120333,  -0.0725717, 
-0.0393023,  -0.0207696,  -0.015532} 
dr0*rho*c''2/ {Omega*al)  =  0.182827  ul*rho*c/Omega=  0.654286 

c*tt/al  =  0.85 

u*rho*c/Omega  =  (1.34624,  1.23551,  1.10859,  0.915125,  0.544341, 

-0.0171408,  -0.393764,  -0.285228,  -0.179471,  -0.10316,  -0.0544994, 
-0.0283824,  -0.0209926} 

strain*rho*c^2 /Omega  =  (-0.241845,  -0.316134,  -0.317649,  -0.249897, 
-0.281069,  -0.259464,  -0.248512,  -0.208936,  -0.152945,  -0.102151, 
-0.0604353,  -0.0356691,  -0.028224} 
dr0*rho*c''2/ (Omega*al)  =  0.197325  ul*rho*c/Omega=  0.677174 

c*tt/al  =  0.9 

u*rho*c/Omega  =  (1.36015,  1.22377,  1.08697,  0.922142,  0.567654, 

0.0837728,  -0.455844,  -0.351223,  -0.239544,  -0.150451,  -0.0860863, 
-0.0494446,  -0.0385668} 

s train* rho*c'' 2 /Omega  =  (-0.256918,  -0.308549,  -0.319547,  -0.267719, 
-0.278596,  -0.262913,  -0.256067,  -0.22669,  -0.183197,  -0.134286, 
-0.0871144,  -0.0572893,  -0.047561} 
dr0*rho*c"2/ (Omega*al)  =  0.211184  ul*rho*c/Omega=  0.703284 

c*tt/al  =  0.95 

u*rho*c/Omega  =  (1.36019,  1.2219,  1.07435,  0.917561,  0.586764,  0.159416, 
-0.507742,  -0.411467,  -0.303403,  -0.206664,  -0.128171,  -0.0809948, 
-0.0661173} 

strain*rho*c^2 /Omega  =  (-0.286447,  -0.295866,  -0.308943,  -0.28983, 
-0.279066,  -0.267091,  -0.255019,  -0.234941,  -0.207834,  -0.165712, 
-0.118435,  -0.0863688,  -0.0747836} 
dr0*rho*c''2/ (Omega*al)  =  0.224143  ui*rho*c/Omega=  0.733475 

c*tt/al  =  1. 

u*rho*c/Omega  =  (1.34356,  1.22935,  1.07396,  0.903457,  0.601409,  0.214984, 
-0.548124,  -0.462743,  -0.366573,  -0.268742,  -0.180898,  -0.125228, 
-0.10638} 

strain*rho*c''2/Omega  =  (-0.320901,  -0.282446,  -0.291901,  -0.309826, 
-0.281019,  -0.273354,  -0.2489,  -0.23612,  -0.224602,  -0.193186, 
-0.152668,  -0.122611,  -0.110265} 
dr0*rho*c"2/ (Omega*al)  =  0.236018  ul*rho*c/Omega=  0.768066 

c*tt/al  =1.05 

u*rho*c/Omega  =  (1.31431,  1.24329,  1.08417,  0.884827,  0.612767,  0.253996, 
-0.383049,  -0.504367,  -0.424746,  -0.332719,  -0.243241,  -0.183366, 
-0.161404} 
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strain*rho*c^2/0mega  =  {-0.348104,  -0.272466,  -0.276778,  -0.32202, 
-0.281378,  -0.281809,  -0.246109,  -0.234177,  -0.233007,  -0.214507, 
-0.187483,  -0.164381,  -0.153053} 
dr0*rho*c''2/ (Omega*al)  =  0.247058  ul*rho*c/Omega=  0.806907 


c*tt/al  =  1.1 

u*rho*c/Omega  =  (1.28199,  1.25959,  1.09936,  0.86796,  0.623299,  0.279507, 
-0.257114,  -0.538117,  -0.474942,  -0.394744,  -0.313019,  -0.255083, 
-0.231854} 

strain*rho*c''2/Omega  =  (-0.358095,  -0.268705,  -0.270472,  -0.323323, 
-0.277993,  -0.291019,  -0.252953,  -0.23334,  -0.235036,  -0.229302, 
-0.220371,  -0.208705,  -0.200613} 
dr0*rho*c''2/ (Omega*al)  =  0.257548  ul*rho*c/Omega=  0.849457 

c*tt/al  =  1.15 

u*rho*c/Omega  =  (1.25812,  1.27401,  1.11261,  0.858416,  0.635252,  0.294992, 
-0.165601,  -0.567473,  -0.516657,  -0.452171,  -0.387183,  -0.338151, 
-0.316395} 

strain*rho*c^2 /Omega  =  (-0.347121,  -0.271834,  -0.275379,  -0.314141, 
-0.271591,  -0.298435,  -0.266084,  -0.236855,  -0.23501,  -0.239227, 
-0.249194,  -0.251635,  -0.248907} 
dr0*rho*c''2/ {Omega*al)  =  0.267283  ul*rho*c/Omega=  0.894874 

c*tt/al  =  1.2 

u*rho*c/Omega  =  (1.25172,  1.28349,  1.11881,  0.859329,  0.649089,  0.304607, 
-0.102485,  -0.596414,  -0.552207,  -0.504322,  -0.46236,  -0.428467, 
-0.41138} 

strain*rho*c^2 /Omega  =  (-0.319097,  -0.280381,  -0.288451,  -0.298047, 
-0.265585,  -0.301381,  -0.280374,  -0.246084,  -0.237979,  -0.247464, 
-0.272723,  -0.28897,  -0.292929} 
dr0*rho*c^2/ (Omega*al)  =  0.276121  ul*rho*c/Omega=  0.9421 

c*tt/al  =  1.25 

u*rho*c/Omega  =  (1.26599,  1.28687,  1.11694,  0.870586,  0.662871,  0.312767, 
-0.0603516,  -0.628159,  -0.585945,  -0.552615,  -0.535529,  -0.520552, 
-0.511028} 

strain*rho*c^2 /Omega  =  (-0.28396,  -0.291352,  -0.302743,  -0.280408, 
-0.264359,  -0.298225,  -0.291068,  -0.260237,  -0.247764,  -0.2576, 
-0.290919,  -0.317165,  -0.327658} 
dr0*rho*c''2/ (Omega*al)  =  0.283971  ul*rho*c/Omega=  0.989952 

c*tt/al  =  1.3 

u*rho*c/Omega  =  (1.29719,  1.28509,  1,11037,  0.889117,  0.672974,  0.323237, 
-0.0315234,  -0.545523,  -0.622659,  -0.599946,  -0.60459,  -0.608452, 
-0.608188} 

strain*rho*c-'2 /Omega  =  (-0.253658,  -0.301255,  -0.310713,  -0.266525, 
-0.270902,  -0.289269,  -0.295808,  -0.278935,  -0.26548,  -0.272266, 
-0.304797,  -0.334182,  -0.349266} 
dr0*rho*c^2/ (Omega*al)  =  0.290966  ul*rho*c/Omega=  1.03208 

c*tt/al  =  1.35 

u*rho*c/Omega  =  (1.33636,  1.28066,  1.1052,  0.910104,  0.675836,  0.338066, 
-0.00979078,  -0.471102,  -0.665801,  -0.649527,  -0.668601,  -0.686849, 
-0.695628} 

strain*rho*c''2/Omega  =  (-0.237511,  -0.307185,  -0.307742,  -0.259956, 
-0.284993,  -0.277206,  -0.295624,  -0.302616,  -0.289123,  -0.29223, 
-0.315846,  -0.340024,  -0.35624} 
drO*rho*c''2/(Omega*al)  =  0.297452  ul*rho*c/Omega=  1.07261 

c*tt/al  =  1.4 

u*rho*c/Omega  =  (1.37284,  1.27681,  1.10727,  0.928633,  0.669898,  0.356684, 
0.00814026,  -0.426323,  -0.716211,  -0.703714,  -0.727575,  -0.752115, 
-0.76757} 

strain*rho*c''2 /Omega  =  (-0.238862,  -0.307618,  -0.294296,  -0.261615, 
-0.302845,  -0.266794,  -0.2944,  -0.325459,  -0.31435,  -0.316183, 
-0.325269,  -0.336785,  -0.35004} 
dr0*rho*c^2/ (Omega*al)  =  0.303339  ul*rho*c/Omega=  1.1115 

c*tt/al  =  1.45 
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u*rho*c/Omega  =  {1.39832,  1.27646,  1.11941,  0.941156,  0.656772,  0.375792, 
0.0222441,  -0.402573,  -0.771835,  -0.762999,  -0.781917,  -0.803012, 
-0.821062} 

strain*rho*c''2 /Omega  =  {-0.254345,  -0.302723,  -0.275778,  -0.269918, 
-0.318511,  -0.263211,  -0.296988,  -0.34254,  -0.336127,  -0.340741, 
-0.333352,  -0.328171,  -0.334951} 
dr0*rho*c''2/ (Omega*al)  =  0.308559  ul*rho*c/Omega=  1.14761 

c*tt/al  =  1.5 

u*rho*c/Omega  =  {1.4095,  1.28138,  1.13991,  0.946273,  0.640975,  0.39034, 
0.0301133,  -0.391146,  -0.828514,  -0.825409,  -0.831749,  -0.840862, 
-0.85672} 

strain*rho*c^2 /Omega  =  {-0.275917,  -0.294148,  -0.260174,  -0.281783, 
-0.326463,  -0.270061,  -0.306775,  -0.351441,  -0.350552,  -0.361157, 
-0.3393,  -0.318591,  -0.317046} 
dr0*rho*c^2/ (Omega*al)  =  0.313077  ul*rho*c/Omega=  1.18005 

c*tt/al  =  1.55 

u*rho*c/Ctaiega  =  {1.40864,  1.29179,  1.16315,  0.944663,  0.628241,  0.395269, 
0.0289221,  -0.38481,  -0.788058,  -0.886642,  -0.876467,  -0.86909, 
-0.87853} 

strain*rho*c''2 /Omega  =  {-0.294637,  -0.28444,  -0.254591,  -0.294001, 
-0.324194,  -0.287882,  -0.323805,  -0.352668,  -0.358195,  -0.372712, 
-0.341683,  -0.312017,  -0.302491} 
dr0*rho*c''2/ (Omega*al)  =  0.317063  ul*rho*c/Omega=  1.20036 

c*tt/al  =  1.6 

u*rho*c/Omega  =  {1.40174,  1.30641,  1.18185,  0.938311,  0.623084,  0.38738, 
0.0175509,  -0.37894,  -0.75965,  -0.941054,  -0.914826,  -0.892221., 
-0.892666} 

strain*rho*c^2 /Omega  =  {-0.304488,  -0.276325,  -0.262226,  -0.3044, 
-0.313616,  -0.313733,  -0.344517,  -0.349007,  -0.361765,  -0.372357, 
-0.339498,  -0.310915,  -0.295707} 
dr0*rho*c''2/ (Omega*al)  =  0.320705  ul*rho*c/Omega=  1.21755 

c*tt/al  =  1.65 

u*rho*c/Omega  =  {1.39539,  1.32281,  1.19004,  0.929464,  0.626716,  0.366643, 
-0.00214865,  -0.371813,  -0.738576,  -0.98329,  -0.945661,  -0.914601, 
-0.905638} 

strain*rho*c^2 /Omega  =  {-0.304613,  -0.272078,  -0.281121,  -0.312371, 
-0.300546,  -0.342,  -0.363091,  -0.344124,  -0.362005,  -0.360028, 
-0.333183,  -0.315593,  -0.298069} 
dr0*rho*c^2/ {Omega*al)  =  0.323957  ul*rho*c/Omega=  1.23126 

c*tt/al  =  1.7 

u*rho*c/Omega  =  {1.39378,  1.33805,  1.18544,  0.919746,  0.636259,  0.336476, 
-0.0253831,  -0.363971,  -0.721675,  -1.01008,  -0.968772,  -0.939186, 
-0.922405} 

strain*rho*c"2 /Omega  =  {-0.299123,  -0.273147,  -0.305186,  -0.318659, 
-0.292407,  -0.366134,  -0.373636,  -0.340982,  -0.359927,  -0.339067, 
-0.324775,  -0.324175,  -0.307695} 
dr0*rho*c^2/ (Omega*al)  =  0.326792  ul*rho*c/Omega=  1.24136 

c*tt/al  =  1.75 

u*rho*c/Omega  =  {1.39698,  1.34924,  1.17012,  0.909771,  0.645721,  0.302882, 
-0.045742,  -0.356965,  -0.706917,  -0.97434,  -0.985484,  -0.966768, 
-0.945095} 

strain*rho*c^2 /Omega  =  {-0.294747,  -0.280042,  -0.326984,  -0.324609, 
-0.295082,  -0.380682,  -0.372553,  -0.340725,  -0.356102,  -0.316372, 
-0.317322,  -0.333261,  -0.320385} 
dr0*rho*c^2/ (Ctoega*al)  =  0.329275  ul*rho*c/Omega=  1.24238 

c*tt/al  =  1.8 

u*rho*c/Omega  =  {1.40138,  1.35402,  1.14934,  0.899324,  0.648384,  0.272769, 
-0.0573476,  -0.352027,  -0.692809,  -0.928424,  -0.998568,  -0.995912, 
-0.972791} 

strain*rho*c"2 /Omega  =  {-0.297511,  -0.292434,  -0.340971,  -0.331273, 
-0.310255,  -0.382895,  -0.360115,  -0.342532,  -0.350412,-0.300277, 
-0.313733,  -0.339184,  -0.331364} 
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dr0*rho*c^2/ (Omega*al)  =  0.331543  ul*rho*c/Omega=  1.2403 


c*tt/al  =  1.85 

u*rho*c/Omega  =  {1.40167,  1.35086,  1.12909,  0.887912,  0.63967,  0.252012, 
-0.0569791,  -0.349258,  -0.677967,  -0.886336,  -1.0115,  -1.0237, 
-1.00239} 

strain*rho*c''2 /Omega  =  {-0.309994,  -0.309335,  -0.345727,  -0.338793, 
-0.334421,  -0.373414,  -0.340557,  -0.344396,  -0.34239,  -0.295352, 
-0.315535,  -0.339413,  -0.337135} 
dr0*rho*c''2/ (Omega*al)  =  0.333575  ul*rho*c/Omega=  1.23642 

c*tt/al  =  1.9 

u*rho*c/Omega  =  {1.39347,  1.33917,  1.11347,  0.875391,  0.619249,  0.243912, 
-0.0452886,  -0.347538,  -0.66111,  -0.853966,  -1.00971,  -1.04696, 
-1.03018} 

strain*rho*c''2 /Omega  =  {-0.33034,  -0.329272,  -0.344227,  -0.346306, 
-0.360071,  -0.355868,  -0.320571,  -0.344219,  -0.331944,  -0.302715, 
-0.322314,  -0.333443,  -0.336694} 
dr0*rho*c''2/ (C)mega*al)  =  0.335372  ul*rho*c/Omega=  1.22858 

c*tt/al  =  1.95 

u*rho*c/Omega  =  {1.37538,  1.31937,  1.10316,  0.862328,  0.591466,  0.248476, 
-0.0265827,  -0.345076,  -0.641472,  -0.835444,  -1.00608,  -1.06355, 
-1.05346} 

strain*rho*c^2/Omega  =  {-0.353348,  -0.350414,  -0.342067,  -0.352386, 
-0.378615,  -0.335612,  -0.306809,  -0.340864,  -0.320063,  -0.319565, 
-0.332113,  -0.322972,  -0.331711} 
dr0*rho*c''2/ {Omega*al)  =  0.337001  ul*rho*c/Omega=  1.2189 

c*tt/al  =  2. 

u*rho*c/Omega  =  {1.34955,  1.29282,  1.09554,  0.849943,  0.563842,  0.262694, 
-0.00722563,  -0.340284,  -0.619411,  -0.831483,  -1.00837,  -1.06822, 
-1.07162} 

strain*rho*c''2 /Omega  =  {-0.372987,  -0.370667,  -0.344565,  -0.35578, 
-0.383853,  -0.318091,  -0.303398,  -0.334722,  -0.309069,  -0.339953, 
-0.340682,  -0.311423,  -0.325643} 
dr0*rho*c''2/ (Omega*al)  =  0.33846  ul*rho*c/Omega=  1.20831 

“  2*05 

u*rho*c/Omega  =  {1.32057,  1.26179,  1.08642,  0.839658,  0.544211,  0.281617, 
0.00670417,  -0.332576,  -0.596792,  -0.838859,  -1.01337,  -1.06424, 
-1.08449} 

strain*rho*c^2 /Omega  =  {-0.385068,  -0.387828,  -0.35407,  -0.356078, 
-0.374536,  -0.307382,  -0.310499,  -0.327564,  -0.302054,  -0.3567, 
-0.34397,  -0.302962,  -0.321983} 
dr0*rho*c^2/ (Omega*al )  =  0.339764  ul*rho*c/Omega=  1.19736 


c*tt/al  =  2.1 

u*rho*c/Omega  =  {1.29346,  1.22929,  1.07247,  0.832466,  0.537661,  0.299819, 
0.0115623,  -0.322709,  -0.576766,  -0.851266,  -1.01682,  -1.05773, 
-1.09362} 

strain*rho*c^2 /Omega  =  {-0.388818,  -0.399822,  -0.368763,  -0.354031, 
-0.354937,  -0.305338,  -0.324481,  -0.321771,  -0.301706,  -0.363833, 
-0.339955,  -0.300707,  -0.322574} 

dr0*rho*c^2/ (Omega*al )  =  0.34092  ul*rho*c/Omega=  1.18682 

c*tt/al  =  2.15 

u*rho*c/Omega  =  {1.27164,  1.19873,  1.05318,  0.828435,  0.544528,  0.312776, 
0.00724153,  -0.312515,  -0.562883,  -0.86128,  -1.01509,  -1.05316, 
-1.10207} 

strain*rho*c^2 /Omega  =  {-0.386698,  -0.405053,  -0.38347,  -0.351335, 
-0.333185,  -0.311495,  -0.339607,  -0.319286,  -0.309042,  -0.358663, 
-0.329633,  -0.305796,  -0.326596} 

dr0*rho*c^2/ (Omega*al)  =  0.341934  ul*rho*c/Omega=  1.17726 

c*tt/al  =2.2 

u*rho*c/Qmega  =  {1.25595,  1.17353,  1.03144,  0.826567,  0.560275,  0.317782, 
-0.00317018,  -0.304156,  -0.557809,  -0.862788,  -1.00673,  -1.054, 
-1.11105} 
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strain*rho*c^2 /Omega  =  {-0.382655,  -0.402782,  -0.392086,  -0.349989, 
-0.318034,  -0.323626,  -0.35048,  -0.320788,  -0.32265,  -0.342768, 
-0.316789,  -0.316972,  -0.33094) 
dr0*rho*c^2/ {C)mega*al)  =  0.342814  ul*rho*c/Omega=  1.16913 

c*tt/al  =  2.25 

u*rho*c/Omega  =  {1.24505,  1.15654,  1.01243,  0.825107,  0.577295,  0.314245, 
-0.0149695,  -0.299253,  -0.562166,  -0.853014,  -0.993207,  -1.06167, 
-1.11948} 

strain*rho*c''2 /Omega  =  {-0.379799,  -0.393421,  -0.390489,  -0.351479, 
-0.315428,  -0.338627,  -0.3542,  -0.325452,  -0.338867,  -0.321513, 
-0.306555,  -0.330952,  -0.331815} 
dr0*rho*c''2/ (Omega*al}  =  0.343566  ul*rho*c/Omega=  1.16282 

c*tt/al  =  2.3 

u*rho*c/Omega  =  {1.23691,  1.14948,  1.0014,  0.822211,  0.587891,  0.303426, 
-0.0239363,  -0.298267,  -0.574011,  -0.833457,  -0.978441,  -1.07512, 
-1.12466} 

strain*rho*c''2 /Omega  =  {-0.378691,  -0.378619,  -0.378549,  -0.356112, 
-0.326257,  -0.353388,  -0.351347,  -0.3314,  -0.352938,  -0.302234, 
-0.303322,  -0.343487,  -0.326813} 
dr0*rho*c^2/ (Omega*al)  =  0.3442  ul*rho*c/Omega=  1.1586 

c*tt/al  =  2.35 

u*rho*c/Omega  =  {1.23051,  1.15256,  1.00125,  0.8167,  0.587218,  0.287823, 
-0.0280455,  -0.300443,  -0.589274,  -0.809359,  -0.967295,  -1.09119, 
-1.12396} 

strain*r}io*c"2 /Omega  =  {-0.377153,  -0.361051,  -0.360339,  -0.362784, 
-0.346198,  -0.365417,  -0.34541,  -0.336632,  -0.360732,  -0.291728, 
-0.308857,  -0.35077,  -0.31648} 
dr0*rho*c^2/ (Qmega*al)  =  0.344728  ul*rho*c/Omega=  1.15661 

c*tt/al  =2.4 

u*rho*c/Omega  =  {1.2268,  1.16436,  1.0112,  0.808624,  0.574998,  0.270469, 
-0.0280858,  -0.304303,  -0.603049,  -0.787907,  -0.963596,  -1.10569, 
-1.11656} 

strain*rho*c^2 /Omega  =  {-0.371662,  -0.34392,  -0.342417,  -0.36928, 
-0.367632,  -0.373117,  -0.340987,  -0.340004,  -0.36033,  -0.293922, 
-0.321486,  -0.35071,  -0.304582} 
dr0*rho*c^2/ (Omega*al)  =  0.345161  ul*rho*c/Omega=  1.15688 

c*tt/al  =  2.45 

u*rho*c/Omega  =  {1.22837,  1.18202,  1.027,  0.7994,  0.55539,  0.254318, 
-0.0269805,  -0.308474,  -0.611293,  -0.775767,  -0.968538,  -1.11482, 
-1.10462} 

strain*rho*c^2 /Omega  =  {-0.359658,  -0.330273,  -0.330982,  -0.372997, 
-0.382777,  -0.375798,  -0.341584,  -0.341763,  -0.352812,  -0.308553, 
-0.336707,  -0.343647,  -0.296819} 
dr0*rho*c''2/ (Omega*al )  =  0.345511  ul*rho*c/Omega=  1.15934 

c*tt/al  =2.5 

u*rho*c/Omega  =  {1.23795,  1.20181,  1.04273,  0.791454,  0.535139,  0.241807, 
-0.0281745,  -0.312425,  -0.612241,  -0.776831,  -0.980159,  -1.11647, 
-1.09307} 

strain*rho*c^2 /Omega  =  {-0.341538,  -0.322307,  -0.329202,  -0.371852, 
-0.386699,  -0.373561,  -0.34799,  -0.343394,  -0.341862,  -0.331324, 
-0.34906,  -0.332262,  -0.298436} 
dr0*rho*c^2/ (Omega*al)  =  0.345792  ul*rho*c/Omega=  1.16379 
DONE 
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